/*
 * TotalEquilibriumNewtonSolver.cpp
 *
 *  Created on: 4 Sep 2011
 *      Author: Allan
 */

#include "TotalEquilibriumNewtonSolver.h"

// C++ includes
#include <cmath>

TotalEquilibriumNewtonSolver::TotalEquilibriumNewtonSolver() :
tolerance(1.0E-6), maxiterations(100)
{}

TotalEquilibriumNewtonSolver::TotalEquilibriumNewtonSolver(const Multiphase& multiphase, const ReactionSystem& eReactions, const vector<string>& jSpecies, const vector<string>& eSpecies) :
tolerance(1.0E-6), maxiterations(100)
{
	Initialiaze(multiphase, eReactions, jSpecies, eSpecies);
}

void TotalEquilibriumNewtonSolver::Initialiaze(const Multiphase& multiphase, const ReactionSystem& eReactions, const vector<string>& jSpecies, const vector<string>& eSpecies)
{
	// Set the number of primary (Nj) and equilibrium (Ne) species
	Nj = jSpecies.size();
	Ne = eSpecies.size();
	
	// Set the indexes of the primary (jIndexes) and equilibrium (eIndexes) species
	jIndexes = multiphase[jSpecies];
	eIndexes = multiphase[eSpecies];
	
	// Initialiaze the total equilibrium problem instance
	problem.Initialiaze(multiphase, eReactions, jSpecies, eSpecies);
	
	// Initialiaze the multiroot solver
	solver.Initialiaze(Nj);
	solver.SetConvergenceOptions(tolerance, maxiterations);
}

void TotalEquilibriumNewtonSolver::SetSolverOptions(double tolerance, unsigned maxiterations)
{
	// Set data members
	this->tolerance = tolerance;
	this->maxiterations = maxiterations;
	
	solver.SetConvergenceOptions(tolerance, maxiterations);
}

void UnderRelaxCHESS(VectorXd& n, const VectorXd& nOld)
{
	const double a =  0.5;
	const double b =  3.0;
	const double c = -0.9;
	
	for(unsigned i = 0; i < unsigned(n.rows()); ++i)
	{
		const double delta = n[i] - nOld[i];
		const double alpha = abs(delta)/nOld[i];
		
		if(alpha > a)
		{
			double r = (delta > 0) ? (b * alpha - a*a)/(b + alpha - 2*a) : c*(alpha - a*a)/(1 + alpha - 2*a);
			
			r *= nOld[i]/delta;
			
			n[i] = nOld[i] + r * delta;
		}
	}
}

const TotalEquilibriumNewtonSolver::Diagnostic TotalEquilibriumNewtonSolver::Solve(const VectorXd& uj, double T, double P, VectorXd& n)
{
	// Set the paramenters of the total equilibrium problem
	problem.SetParameters(uj, T, P, n);
	
	// The vector of moles of primary species (nj)
	VectorXd nj(Nj);
	
	// Initialise the vector (nj) with the values from the initial guess (n)
	for(unsigned j = 0; j < Nj; ++j) nj[j] = n[jIndexes[j]];
	
	// Set the diagnostic instance
	Diagnostic diagnostic = {0, false};
	
	// The solution vector (nj) at the previous iteration (used for application of under relaxation to the solution)
	VectorXd njOld = nj;
	
	// Solve the equilibrium problem
	while(solver.Iterate(problem, nj) && diagnostic.iterations < maxiterations)
	{
		// Apply the under relaxation to the vector (nj) to avoid negative values of compositions
		UnderRelaxCHESS(nj, njOld);
		
		// Update the vector (njOld)
		njOld = nj;
		
		// Update (iteration)
		++diagnostic.iterations;
	}
	
	// Update the composition of the equilibrium species for the last time
	problem.UpdateEquilibriumComposition(nj);
	
	// Transfer the solution (nj) and (ne) to the vector (n) (use the equilibrium problem instance for this)
	n = problem.n;
	
	// Check convergence of the previous iterations
	diagnostic.converged = diagnostic.iterations < maxiterations;
	
	return diagnostic;
}
